/*
 * Interpolate.h
 *
 *  Created on: May 30, 2010
 *      Author: snirgaz
 */

#ifndef Interpolate_H_
#define Interpolate_H_

#include "HelpLibs/def.h"
#include "ParticlePos.h"
#include "SimParams.h"
#include "WorldLines.h"
#include "HelpLibs/hdf5RW.h"
#include <vector>

using namespace std;

enum O_P {
	ORIGINAL, PROPOSED
};
#ifdef POT_AZIZ
void calcGreenForce(double dr, double drTh, double rZero, double &v, double &f);
void calcGreen(double dr, double drTh, double rZero, double &v);
#endif
//class Interpolate: virtual public QmcParams {
//public:
//	Interpolate();
//	~Interpolate();
//	void calcEnergy();
//	void calcDist();
//	void calcVF();
//	int interpIndex(BeadPos const &Pos, int const &res);
//	void calcSingleGreen(BeadPos const &Pos, double &v, ForceData &f);
//	void calcSingleGreen(BeadPos const &Pos, double &v);
//#ifdef EXT_POT
//	void calcEnergyExtPot(BeadPos Pos, double &v, ForceData &f);
//	void calcEnergyExtPot(BeadPos Pos, double &v);
//#endif
//
//protected:
//
//#ifdef EXT_POT
//#ifndef LATTICE
//	__declspec(align(16)) double *vgridExtPot_;
//	std::vector<double *> fgridExtPot_;
//#endif
//	std::vector<ForceData> extPotFOriginal_, extPotFProposed_;
//	ForceData extPotJumpFOriginal_, extPotJumpFProposed_;
//#endif EXT_POT
//	vector<double> resV_, resF_;
//	O_P OP_;
//#ifdef POT_LOG
//	__declspec(align(16)) double *LogArg_;
//#endif
//#ifdef EXT_POT
//	double dResExtPot_, dxExtPot_;
//#endif
//#ifdef POT_AZIZ
//	vector<double> drSqr_;
//#endif
//};
//#ifdef POT_LOG
//template<int N >
//void calcGreen(double *restrict x, double * restrict y, double * restrict v,
//		double *restrict logArg, int size) {
//	double xArg, yArg, invArg;
//	double trigArg[size], hypArg[size], c[N][size], ch[N][size], s[N][size],
//	sh[N][size];
//	//#pragma parallel
//	//#pragma loop count min(16)
//#ifdef ROT
//	double const t[3] = {-0.008629446001672822, 0.00001877888595624988,
//		-5.4250183199002677E-8};
//	double const cost = 0.5, sint = 0.8660254037844386;
//	trigArg[0:size]=2 * PI * (y[0:size] + x[0:size] * cost);
//	hypArg[0:size]=2 * PI * sint * x[0:size];
//#else
//	double const t[3] = {-0.008629446001672822, 0.00001877888595624988,
//		-5.4250183199002677E-8};
//	double const cost = 0.5, sint = 0.8660254037844386;
//	trigArg[0:size]=2 * PI * (y[0:size]);
//	hypArg[0:size]=2 * PI * x[0:size];
//#endif
//	ippsCos_64f_A26(trigArg, c[0], size);
//	ippsCosh_64f_A26(hypArg, ch[0], size);
//	if (N >= 2) {
//		c[1][:] = 2 * (c[0][:] * c[0][:]) - 1;
//		ch[1][:] = 2 * (ch[0][:] * ch[0][:]) - 1;
//	}
//	if (N >= 3) {
//		c[2][:] = 2 * c[0][:] * c[1][:] - c[0][:];
//		ch[2][:] = 2 * ch[0][:] * ch[1][:] - ch[0][:];
//	}
//	for (int i = 0; i < size; i++) {
//#ifdef ROT
//		double cpv = PI * sint * pow(x[i], 2);
//#else
//		double cpv = PI * pow(x[i], 2);
//#endif
//		v[i] = cpv;
//		for (int n = 0; n < N; n++)
//		v[i] += t[n] * c[n][i] * ch[n][i];
//		logArg[i] = ch[0][i] - c[0][i];
//	}
//}
//
//template<int N>
//void calcGreenForce(double *restrict x, double *restrict y, double *restrict v,
//		double *restrict logArg, double *restrict fx, double *restrict fy,
//		int size) {
//	double xArg, yArg, invArg;
//	double trigArg[size], hypArg[size], c[N][size], ch[N][size], s[N][size],
//	sh[N][size];
////#pragma parallel
////#pragma loop count min(16)
//#ifdef ROT
//	double const t[3] = {-0.008629446001672822, 0.00001877888595624988,
//		-5.4250183199002677E-8};
//	double const cost = 0.5, sint = 0.8660254037844386;
//	trigArg[0:size]=2 * PI * (y[0:size] + x[0:size] * cost);
//	hypArg[0:size]=2 * PI * sint * x[0:size];
//#else
//	double const t[3] = {0.003741873197321289, 3.487354517808118E-6,
//		4.341608118994279E-9};
//	trigArg[0:size]=2 * PI * (y[0:size] );
//	hypArg[0:size]=2 * PI * x[0:size];
//#endif
//	ippsSinCos_64f_A26(trigArg, s[0], c[0], size);
//	ippsCosh_64f_A26(hypArg, ch[0], size);
//	ippsSinh_64f_A26(hypArg, sh[0], size);
//	if (N >= 2) {
//		c [1][:] = 2 * (c[0][:] * c[0][:]) - 1;
//		ch[1][:] = 2 * (ch[0][:] * ch[0][:]) - 1;
//		s[1][:] = 2 * c[0][:] * s[0][:];
//		sh[1][:] = 2 * ch[0][:] * sh[0][:];
//	}
//	if (N >= 3) {
//		c [2][:] = 2 * c[0][:] * c[1][:] - c[0][:];
//		ch[2][:] = 2 * ch[0][:] * ch[1][:] - ch[0][:];
//		s[2][:] = c[0][:] * s[1][:] + c[1][:] * s[0][:];
//		sh[2][:] = ch[1][:] * sh[0][:] + ch[0][:] * sh[1][:];
//	}
//	for (int i = 0; i < size; i++) {
//		double trigArg, hypArg, cpv, cpf;
//		double c_sh, ch_s;
//
//#ifdef ROT
//		cpv = PI * sint * pow(x[i], 2);
//#else
//		cpv = PI * pow(x[i], 2);
//#endif
//
//		v[i] = cpv;
//		logArg[i] = ch[0][i] - c[0][i];
//		invArg = 1 / logArg[i];
//#ifdef ROT
//		fx[i] = PI
//		* (sqrt(3) * x[i]
//				- 0.5 * (s[0][i] + sqrt(3) * sh[0][i]) * invArg);
//		fy[i] = PI * (-x[i] + 0.5 * (sh[0][i] - sqrt(3) * s[0][i]) * invArg);
//#else
//		fx[i] = 2*PI*(x[i] -0.5 * sh[0][i] * invArg);
//		fy[i] = PI* s[0][i] * invArg;
//#endif
//		for (int n = 0; n < N; n++) {
//			v[i] += t[n] * c[n][i] * ch[n][i];
//			ch_s = ch[n][i] * s[n][i];
//			c_sh = c[n][i] * sh[n][i];
//#ifdef ROT
//			double dfx = PI * t[n] * double(n + 1) * (-ch_s + sqrt(3) * c_sh);
//			double dfy = -PI * t[n] * double(n + 1) * (sqrt(3) * ch_s + c_sh);
//#else
//			double dfx = 2*PI * t[n] * double(n + 1) * c_sh;
//			double dfy = -2*PI * t[n] * double(n + 1) * ch_s;
//#endif
//			fx[i] += dfx;
//			fy[i] += dfy;
//		}
//	}
//}
//#endif
#endif /* Interpolate_H_ */
